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ABSTRACT 



The partial differential equation for heat diffusion is numerically 
integrated by the Runge-Kutta-Glll method. A solution is obtained for 
the diurnal temperature variation with a bounded coefficient of eddy 
diffusivity which varies periodically with time and nonlinearly 
with height. The surface wave is represented by the sum of a 
diurnal and a semidiurnal harmonic wave. The results may be interpreted 
to apply over a fairly broad range of diffusivity values and height. 

With appropriate choices of the various parameters, reasonably good 
agreement is obtained between theoretical and observational values of 
amplitude reduction and phase lag. 

The author wishes to express his appreciation for the assistance 
and encouragement of Professors G. J. Haltiner and E. J. Stewart of 
the U. S. Naval Postgraduate School in this investigation. 
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1 . Int roduc 1 1 or. . 

Of the toany problems f&cin ■- meteorologists, jue that mac received 
ranch attention is t?ict of* heat diffusion in the atmosphere. Much of this 
effort has been directed toward the diurnal temperature wave, and in 
every attempt, the concept of eddy diffusivity has been retained. The 
classical theory of daily temperature variations has been discusser 1 by 
Sutton (7), while the more recent contributions have been s acinar l zed 
by Staley (6). 

The Taylor heat -diffusion equation, which states that the 
turbulent flux of heat is proportional to the gradient of potential 
temperature, may be written: 




Mere t represents the- time ; z, the height; and @ (z,t), the potential 
temperature deviation from a mean value. Generally, the coefficient of 
eddy diffusivity, T, is a function of both height and time in any 
part icu'lar location . 

In many of the past studies, K has been assume 1 to vary, (a) 

linearly with height, (b ) as a power of height, and (c) as a bounded 

exponential function of height. In each of these cases, the resumed 

functional form of X was assumed to hold throughout the atmosphere. 

Recently, However, do Vries (l) has studied: 

the case of a- atmospnere consisting of several layers, 
in each of which X is expressed as a different function 
of height. 

For ar example, he proposed three layers; a la.mrer sublayer, a 
turbulent bcund v *v la/^r, and a T ’freo M atmospheric layer. In general, 
each of these layers, -ml ‘ n particular, the turbulent bo Lindary 
lever, can be farther subdivided into layers. ;Ie also gave 
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analytical solutions for several particular cases; (a) K = constant, 
(b) K = o(. ( ^ , and (c) K = ajjL - b exp(-cz)] . Case (a) yields a 
solution In terras of a rather complicated exponential function, (b), 
a solution In terras of Bessel functions of the first kind, and (c), 
a solution in terras of a hypergeoraetric function. However, no 
numerical values were given, and the results are not suitable for 
practical purposes. 

It is the purpose of this investigation to present a numerical 
solution of 2qn. (l), with appropriate boundary conditions, for a 
rarticular case of an atmosphere composed of an arbitrary number of 
layers . 
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2. The Heat Diffusivity Coefficient. 



From physical considerations, K, the heat diffusivity, may be 
expected to increase in the lover layers, but should not continue to 
increase indefinitely with height. Hence, a bounded diffusivity 
coefficient, which increases with height from the surface to some fixed 
level and then decreases with height to some residual value at higher 
levels, would seem more suitable than either the case of K increasing 
linearly with height or of K increasing as a simple power of height. 
Upon examination of Fig. (l), which is based on observations taken 
at the micrometeorological tower at the University of Washington (2), 
(5) and Mildner’s pilot balloon observations taken near Leipzig (5), 
it is seen that no single functional form of K will suffice. The 
values of K from the University of Washington data were obtained 
by methods based on energy continuity at the earth’s surface, 
while those from the Leipzig data were actually values of eddy 
viscosity. These values of eddy viscosity were used as an 
approximation to eddy diffusivity, since they showed the desired 
variation with height, and since values of eddy diffusivity at 
sufficiently great heights are difficult to obtain. Upon 
dividing the atmosphere into four layers, and then fitting functions 
of height to each layer, the following form of K resulted: {see Fig. 



The quantities ,a^,aj,a^,a^.,a^, and h are constants, appropriate 




values of which will be assigned later. The discontinuity which 
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occurs at z = 228 meters is resolved "by defining j £ at this level as zero. 

Haltiner (4) suggests that K should have a diurnal variation, and that 

it appears reasonable to expect that K would more or less follow 
the surface temperature wave, reaching maximum and minimum values 
near the time when the temperature is maximum and minimum, 
respectively. 

The eddy dlffusivity, K, can now be expressed as a function of height times 
a function of time, i. e., K = f(z)g(t), where g(t) is defined: 

g(t) = 1 + b (sin Cot + c sin 2 iOt) (3) 

Here b and c are constants, and CO *= 2TT/86,400 sec * . 

Appropriate boundry conditions associated with Eqn. (l) will now be 
chosen. These are frequently taken to be: 

O - 0 for z =oO y and all t; (4) 

& = C9 0 (t ) for z » 0; (5) 

where 0 o {t) is defined by: 

<%(t) » cU sinCot + c sin 2wt) (6) • 

Here d is a constant and c and OO are the same as in Eqn. (3). 
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J. Transformation of Coordinates. 

Tn order to Diace the basic equations into a form which will give 
a broader interpretation of the results, and to scale these equations 
for computational purposes, let: 

z = lOOq a , and t = 3oOC't . (7 ) 

Here q is a constant, while •'t and are the new variables. When t is 
in seconds, the units of <- are hoiirs, and when z is in centimeters, the 
units of cx are in "q-meters"; i. e., in 1-iieter units when q = 1, 2- 
meter units when q = 2, etc. With this transformation, (7), Fqs. (l) 
and (2) become: 



o) <9 - .<?■ ?A j j u 1 

dt ' r d<s L n 



h - ) 



r 



3 (t) 



cl, cr 

-+ a 



3 



0/»i T-% (or -h)l 



(9) 



Q 
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Here ^ must now be taken as 277 /2h, and; 



a; = a, ( lOOq ) s , 


a \ = b-j ( lOOq ) 


aj. = a 5 ( lOOq ) , 


h’ = h/lOOq. 



The boundary conditions, Eqs. (^) and (5), remain identical in form. 

Both observation and theory indicate that the most pronounced 
variations of temperature occur near the surface of the earth. 
Therefore, a small grid distance is desirable at low levels, while 
at higher levels, a small grid distance is not necessarily needed. 
This suggests that the following transformation of the vertical 
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coordinate may be useful: 



s =• In O' . 



(10) 



With this transformation, the lover boundary condition, Eqn. (5), may 
apply at a height of q meters corresponding to ^ - 1, however, to reduce 
surface effects, it is convenient to select <5 = 2. Actually the vertical 
cooi'd inate system is arbitrary, and can be chosen to have its base at any 
level. By utilising Eqn. (10), "qs. (8) and (9) now become: 








Qj £ X/*[~ % " j J 




( 12 ) 



Here a 7 is taken as I0 e , and the other constants, a’ , ai, , a^,, a y , 
and a^, must be scaled by a factor of 10 -6 . 
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U. Finite Difference Equations. 



In order to obtain a numerical solution, the derivatives of (Q 

with respect to s and -fc in Eqn. (ll) must be replaced by appropriate 

finite difference forms. The problem is then reduced to that of 

solving a system of linear algebraic equations in the values of (9 

J z 6>- 

over a grid of points covering the desired range. Expressing g 

\ g 

and j-g in Eqn. (ll) as finite difference ratios, ve obtain: 



dt 



<9.34 <3 r 

e I (K) 



£ -Ziff; 



J? 




Here the subscript i designates the i ’th level of the vertical grid 
at height i jp , where $ is the vertical distance between grid points in 
units of s, that is, a pure number, since s is dimensionless. 

The partial differential equation, Eqn. (l), has now been reduced 
to a system of ordinary differential equations which will be 
integrated numerically by the Runge-Kutta-Gill method. 
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5. The Const ante. 

Kext, appropriate values or./', n, anoxic, as well as the constants 
a* , nj, etc., will now be chosen. ince the vertical coordinate is 
s = lr. <5 , $ - /] s, a constant . '-’ith a choice of J( - In 2, the i 'th 
eauation of the system defined by Fqn. 'll) applies at the height 
2 q- meters. Therefore, the successive values of i, beginning with 
the lower boundary condition, aoply to the levels, 2, 4, 8, In, 12,64, 
etc. q-me + ers. To reduce the integration time, yet still give a fairly 
detailed nicture of the vertical structure, the upper boundary condition, 
'•’qn . f4), Is chosen to apply at the height of 1024 q-meterc . For this 
choice, & ~0; an 0 the system, Eqn (id), to be solved consists of 

eight simultaneous ordinary differential equations, uojether with 
Fqn. (6) for <5^, . A numerical solution for this system was then 
obtained by the Runge-Kutta-Oill method on a Kational Cash Register 
102A electronic computer. The choice of an appropriate time interval 
depended intimately upon the size of the eddy diffusivity. As the 
latter increased , A b had to be decreased in order to maintain 
computational stability and reasonable accuracy. 

The following muierical values of the other constants were 
chosen as being representative: 



a' - 0.02 (sec ' ), 

- 52,500 (cm 0 /sec ) , 
a^V = 0.0085, 
b =■ 0.3, 

6 = 0.05 (°C), 



a L = 205 0 (cm/sec), 

= 415,000 (cm 2 /sec), 
= 53,000 (cn 2 /sec), 
c = 0.5, 
u* - 223. 



( 14 ) 



r he choice of a) , , a.,, a, /} and a, ; gives a >28- fold increase 

of K with height from the surface to 223 a-meters, and a 14-fold 
decrease of K with height above 228 q-meterc. Also, the cnoice of 
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b and c gives rise to a 21-fold diurnal variation of K with time. The 
value of d was chosen because it keeps the magnitude of$ o <0.1. 

With the value of 10 s , the integration interval was 

prohibitively short for the computer used; therefore a smaller 
value of ay = 6.25 x 10 4 was used. This permitted a time interval 
of 1/3O hour. 

The value of a 7 = 6.25 x 10 4 gives the following range of the 
coefficient of eddy diffusivity, Eqn . (12), in units of cm 2 /sec: 



at surface. 


K mx - 0.15 x 10 2 q 2 , 


^/7t/ A/ 


= 0.714 q : 


at 20 q-m, 


= °- 15 x 1C)3< 1 ? > 




* 7.14 q 2 


at 228 q-n, 


K m/ ■ 0 - 4 9 2 * 10 


^ miKi 


= 230 q 2 , 


at 460 q-n, 


K /W/J/ = °’35 x 10 3 q 2 , 




= 17 q 2 . 



For q < 1, the values of K are somewhat small for typical atmospheric 
conditions, esnecially in the lowest layers. This is reflected in 
the results by a large amplitude reduction and phase lag with 
increasing height. However, such values of K are not unrealistic 
for winter. 
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r y, Results. 

The results of the integration are given hy Fir. ( 3 ) and Table (l), 

Fig. (‘5) being the graphical interpretation of Table (l). Ir. Fi^. (3), 
amplitude reduction and phase lag, in hours, are plotted against height, 
in q-metors. As car be seen fr on FI'-. (3), the amplitude reduction and 
chase lag are quite large, especially when compared with observed values. 

(Fee Pig. (h) which is based on observations taken at the microoieteorological 
tower at the University of Washington ) . The large amplitude reduction 
and phase lag from the theoretical results. Pig. ( 3 ), are due largely 
to the reduction of a from 10 6 to '.25 x 10 4 . However, by proper 
selection of the value of q, fair agreement between theoretical and 
observed values can be obtained. Fig. ( 3 ) also shows the amplitude 
reduction and phase lag for the minimum at three levels. Fro n these 
three levels, it can be seen that the amplitude reduction is less 
for the maximum than for the minimum, the latter corresponding to the 
lower nocturnal diffusivlty . This difference is appreciable in that at 
the eight q-meter level, the maximum is about ^ 3 °/o of the surface value 
o- Q while the minimum is only 19 u /o of the surface minimum. At lower 
levels, however, this difference is much smaller. The phase lag increases 
with heigh' 1 for both the maximum and the minimum, however, the lag 
of the maximum is greater than that of the minimum at any particular 
level. Tli i s does not seem consistent with greater amplitude reduction of 
minima compared to maxima. Possibly, these differences would vary 
in succeeding maxima and minima, as greater time elapses from the 
initial conditions. 

By a proper choice of q, some comparisons between observed and 
theoretical results may be made. Table (2) shows amplitude reduction 
and phase lag from observations taker., at Ieafield, along with 



12 



MEAN RELATIVE AMPLITUDE AND PHASE LAG 
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theoretical values for q = 10 and q = 20. For both values of q, there 

is fair agreement with observed values, however, both the theoretical 

amplitude end phase lag increase with height much more rapidly than do 

the observed values. Table (3) compares theory with observed data at 

the Eiffel Tower. Here q is taken to be 60 and 33, the value of 60 

comparing to the summer data and the 33, to the winter data. Table (4) 

shows some observations taken at the University of Washington. By the 

choice of q = 8, the amplit\ide compares fairly wellj but where the observed 

phase lags are the order of minutes, those from theory are the order 

of hours. In Table (5), the data from Forton i3 compared to theory. For 

# 

the choice of q = 2.5, the comparison is fair, but once again, shows 
the theoretical phase lag and amplitude reduction much greater than 
the observed. Table (6) compares the amplitude reduction in summer 
and winter at Ieafield with theory. Here q = 5 was chosen for winter, 
and a = 12, for summer. For the winter case the theoretical and 
observed values compare very well, but for summer, the theoretical 
amplitude reduction is significantly greater than the observed value. 
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7 . Conclusions . 



A numerical solution has teen presented for the diurnal temperature 
variation with a coefficient of eddy diffusivity which is a function 
of height and time. By suitable selection of several parameters, 
reasonably good agreement has been obtained between theory and 
observation . 

Improvement in the results may possibly be achieved by improving 
the functional form of the diffusivity coefficient, especially in the 
lowest layers. It would also be desirable to include such parameters as 
surface roughness, stability, wind velocity, etc. With a diffusion 
coefficient in terms of these parameters, as well as height and time, 
a solution for the diurnal temperature wave would have a much broader 
application. 
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